風紋の数理モデル¶


砂丘とかでみられる波波のやつ、風紋の数理モデルがあるらしいので簡単に実装してみた。

風紋モデルの運動には3つのステップがあるようで、saltation, suspension, creep という感じらしい。

saltationは日本語で跳躍とか跳動とか訳されるようで、砂粒が風で飛ばされてしばらくして着陸するみたいな運動。今回のモデルではランダムな点で決まった回数だけ跳躍させた。

suspension(浮遊)は黄砂みたいなフワフワと長く飛ばされる現象で、扱いが面倒なのでモデルには入れない。

creepはそのまんまクリープで、風やら重力やらで表面を転がるのを表している。このモデルでは高さの拡散効果みたいなもんだと思う。

・saltation¶

二次元の高度の場を $h(x,y)$とする。

風は$X$正の方向に吹くとして、砂の飛距離を、$L=L_0+bh(x,y)$ で表す。$L_0$は風の効果、$bh$は高さの効果(高いところの砂粒のほうが飛ぶでしょということ)。

あとは簡単で、飛ばされた砂の量(高さ)を $q$とすると、saltation が起こったセル(砂が運ばれる前と後のセル)では

$$ h'(x,y)=h(x,y)-q, \quad h'(x+L(h(x,y)),y)=h(x+L(h(x,y)),y)+q $$

と高度が変化する。$h'$は、次の時間ステップに行く前の中間ステップ。

・creep¶

拡散方程式を有限差分法で解くみたいなもん。省略。


上のふたつを交互に発展させて計算するぞ!初期値は適当にランダムノイズを乗せとく。

境界条件は周期境界、計算領域は100$\times$100で、時間発展は1000ステップにしておいた。

In [1]:
using Plots
In [2]:
function main()
    dh=0.3 
    nx=100
    ny=100
    q =0.1
    b =1.0
    L0=2.0
    D =0.1

    h=dh*rand(ny,nx).+2.70 # initial condition
    hnew=copy(h)
    N=1000 # timestep
    freq=10
    
    Hs = Array{Array{Float64,2}}(undef, N÷freq+1)
    
    num=1
    Hs[num]=copy(h)

    for itr=1:N

        # saltation
        for k=1:nx*ny
            ir=rand(1:nx)
            jr=rand(1:ny)

            if h[jr,ir]>0
                L=Int(floor(L0+b*h[jr,ir]+0.5))
                ib=ir+L
                h[jr,ir]-=q
                h[jr,ifelse(ib>nx, ib%nx, ib)]+=q
            end
        end
        
        # creep
        for i=1:nx
            for j=1:ny
                hN=h[ifelse(j==ny,1,j+1),i]
                hS=h[ifelse(j==1,ny,j-1),i]
                hE=h[j,ifelse(i==nx,1,i+1)]
                hW=h[j,ifelse(i==1,nx,i-1)]

                hNE=h[ifelse(j==ny,1,j+1),ifelse(i==nx,1,i+1)]
                hSE=h[ifelse(j==1,ny,j-1),ifelse(i==nx,1,i+1)]
                hSW=h[ifelse(j==1,ny,j-1),ifelse(i==nx,1,i+1)]
                hNW=h[ifelse(j==ny,1,j+1),ifelse(i==1,nx,i-1)]

                hnew[j,i]=(1-D)*h[j,i]+D*((hN+hS+hE+hW)/6+(hNE+hSE+hSW+hNW)/12)
            end
        end
        
#         for i=1:nx
#             for j=1:ny
#                 h[j,i]=hnew[j,i]
#             end
#         end
        @.(h=copy(hnew))
        
        if itr%freq==0
            num+=1
            Hs[num] = copy(h)
            print(itr, " ")
        end
    end

    return Hs
end
@time H=main();
10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 940 950 960 970 980 990 1000   0.406244 seconds (33.68 k allocations: 9.676 MiB, 5.38% compilation time)
In [3]:
N=length(H)
anim = @animate for t=1:N
    h=heatmap(H[t],clim=(0,10),c=:sandyterrain,title="height")
    plot(h, size=(500,450))
end
filename="fuumon.gif"
gif(anim, filename, fps=10)
[ Info: Saved animation to C:\****\****\****\****\fuumon.gif
Out[3]:

うおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおお

なんか波ができとりますわ。きれい。

In [4]:
N=length(H)
anim = @animate for t=1:N
    h=surface(H[t], zlims=(0,30),legend=:none,c=:summer)
    plot(h, size=(800,800))
end
filename="fuumon3D.gif"
gif(anim, filename, fps=10)
[ Info: Saved animation to C:\****\****\****\****\fuumon3D.gif
Out[4]:

三次元にするとアスペクト比がおかしいのでギザギザすぎた。

$L$や$q$の扱いをちょっと変えるとバルハンとかの再現もできそうなので、いつかやるかもしれん。

参考¶

西森拓, and 大内則幸. "飛砂による地形の動力学: 風紋と砂丘." 物性研究 61.1 (1993): 32-43.